clear all
%Constants (all MKS, except energy which is in eV)
hbar=1.055e-34;m=9.110e-31;epsil=8.854e-12;q=1.602e-19;

%Lattice
L=5e-10;Np=4;a=(L/Np);% *1 for Fig.1.3.6 and *2 for Fig.1.3.7
R=a*[1:1:Np];
t0=(hbar^2)/(2*m*(a^2))/q
%Quantum numbers
n=1;l=0;% for 1s, n=1 and for 2s, n=2
%Hamiltonian,H = Kinetic,K + Potential,U
K=(2*t0*diag(ones(1,Np)))-(t0*diag(ones(1,Np-1),1))-(t0*diag(ones(1,Np-1),-1));
U=((-q/(4*pi*epsil)./R)+(l*(l+1)*hbar*hbar/(2*m*q))./(R.*R));U=diag(U);
H=K+U;
save Hmat_hydrogen.dat H;
[V,D]=eig(K+U);D=diag(D);[DD,ind]=sort(D);
E=D(ind(1))
psi=V(:,ind(1));
P=psi.*conj(psi);
h=plot(R,P,'b');
xlabel(' Lattice site # --->');% Part (b)
ylabel(' Probability ---> ');% Part (b)
grid on;
save eigen_hydrogen1d.dat D;
Z = V.*V;
save probab_hydrogen1d.dat Z;
save eigen_hydrogen1d.dat V;